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Integration  between  COMSOL  Multiphysics™  and  MATLAB™  offers  a  useful  option  for  the  self-automated 
geometry  optimization  in  proton-exchange  membrane  fuel  cells  (PEMFCS).  It  overcomes  the  difficulties 
of  automatically  re-generating  high-quality  computational  meshes  and  subsequently  running  the  sim¬ 
ulations  to  evaluate  the  objective  function  values  using  commercial  software  in  computational  fuel  cell 
dynamics-based  designs.  Geometry  optimization  studies  of  an  air-breathing  PEMFC  searching  for  the  opti¬ 
mum  channel  ratio  at  the  anode  and  the  optimum  open  ratio  at  the  cathode,  are  undertaken.  A  sequential 
quadratic  programming  method  is  selected  to  deal  with  the  constrained  design  problems,  while  the  objec¬ 
tive  functions  are  evaluated  by  running  the  three-dimensional  simulation  script  of  COMSOL™  under  the 
MATLAB™  environment.  Simulation  results  show  that  for  the  air-breathing  PEM  fuel  cell  operated  at  353  K 
and  one  standard  atmosphere  pressure,  when  the  anode  channel  ratio  is  fixed  at  10%,  the  optimum  cath¬ 
ode  open  ratios  are  very  similar  for  the  cell  operated  at  voltages  of  0.7  and  0.4  V,  namely,  49.8%  for  0.7  V 
and  49.5%  for  0.4  V.  When  the  cathode  open  ratio  is  set  at  80%  with  a  cell  voltage  of  0.7  V,  the  optimum 
anode  channel  ratio  is  found  to  be  34.7%. 

©  2008  Elsevier  B.V.  All  rights  reserved. 


1.  Introduction 

Recently,  air-breathing  proton-exchange  membrane  fuel  cells 
(PEMFCs)  have  been  identified  as  one  of  the  most  promising  alter¬ 
natives  for  traditional  power  sources  in  vehicular  and  portable 
devices.  The  performance  of  a  PEMFC  depends  on  a  variety  of  struc¬ 
tural  and  functional  parameters,  such  as  the  geometry  of  the  flow 
paths  for  the  fuels,  along  with  the  operating  conditions,  thickness  of 
gas  distribution  layers  (GDL),  as  well  as  the  performance  of  the  cata¬ 
lyst.  Studies  indicate  that  the  performance  of  air-breathing  PEMFCs 
is  seriously  affected  by  oxygen  mass  limitation  at  the  cathode  side. 
To  improve  the  performance  of  air-breathing  PEMFCs,  it  is  neces¬ 
sary  to  understand  the  relationship  between  the  oxygen/hydrogen 
mass  fraction  distributions  and  the  geometric  designs  of  the  cell. 
Numerical  optimization  thus  becomes  an  important  tool  for  under¬ 
standing  the  physical  phenomena  related  to  various  geometric 
configurations  that  affect  fuel  cell  performance. 

Unlike  applications  in  structures  and  control  systems,  sensi¬ 
tivity  information  cannot  be  easily  extracted  from  computational 
fluid  dynamics  (CFD)  codes.  Hence,  the  characteristics  of  the  entire 
fuel  cell,  including  the  electrical  field,  the  flow-field  and  chemi¬ 
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cal  reactions,  must  be  computed  in  a  highly  coupled  manner  to 
assess  the  objective  function  values  during  the  fuel  cell  optimiza¬ 
tion  design  cycles.  Consequently  the  efficiency  and  accuracy  of  the 
simulation  are  very  important.  CFD-based  modeling  of  PEMFCs  has 
been  developed  during  the  past  decades  [1,2].  Besides  these  in- 
house  developed  CFD  codes,  there  is  commercial  software,  such 
as  Fluent™,  Star-CD™  and  COMSOL™  that  can  simulate  fuel  cell 
performance  with  reasonable  accuracy. 

The  recent  development  of  CFD  techniques  has  now  made  it 
possible  to  carry  out  CFD-based  optimization  of  the  design  of  fuel 
cells.  There  have  been  several  studies  of  the  effects  of  functional 
and  structural  parameters  on  the  performance  of  the  PEM  fuel  cells. 
Wang  et  al.  [3]  investigated  the  influence  of  electrode  composition 
on  cathode  performance  by  carrying  out  a  series  of  simulations. 
Sinha  et  al.  [4]  found  that  there  exists  an  optimum  value  of  the  dif¬ 
fusive  length  of  diffusion  media  that  provides  peak  performance 
at  a  given  operating  temperature.  Hontanon  et  al.  [5]  presented  a 
study  of  the  performance  of  the  flow  distribution  system  in  terms 
of  the  fuel  consumption  at  the  anode.  Na  and  Gou  [6]  reported  an 
optimization  study  of  the  efficiency  and  cost  of  a  fuel  cell  system 
under  various  operating  conditions.  In  addition  to  examinations  of 
the  effects  of  functional  parameters,  optimization  of  the  structural 
parameters  to  maximize  the  performance  of  PEMFCs  has  also  been 
performed.  Sinha  et  al.  [7]  compared  the  performance  of  PEM  fuel 
cells  with  parallel  and  serpentine  flow-fields  operated  at  elevated 
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Nomenclature 

c  concentration  of  species  (mol  m-3 ) 

D  diffusion  coefficient  (m2  s-1 ) 

Eeq  equilibrium  voltage  (V) 

F  body  force  (N) 

F  Faraday’s  constant  (96,487  A  s  mol-1 ) 

J  transfer  current  density  (A  nrr 2 ) 

jo  exchange  current  density  (A  m-2 ) 

M  molecular  mass  (kg  mol-1 ) 

N  mass  flux  (kg  m-2  s-1) 

p  pressure  (Pa) 

Q  source  term  of  mass  balance  equations  (kg  m-3  s-1 ) 

R  gas  constant  (J  mol-1  K-1 ) 

Ragg  agglomerate  radius  (m) 

s  specific  area  (m-1 ) 

S  current  source  term  (A  m-3 ) 

u  velocity  vector  (ms-1) 

x  mole  fraction 

x  design  variable  vector 

Greek  letters 

r)  dynamic  viscosity  ( N  s  m-2 ) 

rjv  overvoltage  (V) 

k  permeability  (m2) 

A  Lagrange  multiplier 

p  mixture  density  (kg  m-3 ) 

a  conductivity  ( S  m- 1 ) 

0  potentials  (V) 

co  mass  fraction 

Subscripts 
a  anode 

c  cathode 

H2  hydrogen 

H20  water 

m  membrane 

02  oxygen 

s  electrode 

Superscripts 
agg  agglomerate 

ref  reference 

T  temperature 


temperature  and  low  relative  humidity.  Wang  et  al.  [8]  analyzed 
the  effect  of  channel  configuration  on  air-breathing  fuel  cell  per¬ 
formance  by  varying  the  widths  of  channels.  Lin  et  al.  [9]  studied 
the  effect  of  the  channel-to-rib  width  ratio  and  the  porosity  of  the 
GDL  and  the  catalyst  layer  with  a  simplified  conjugate-gradient 
method  based  on  a  two-dimensional  simulation  model.  Hsieh  and 
Chu  [10]  examined  the  effects  of  channel  and  rib  widths  on  cell 
performance  in  terms  of  voltage-current/pressure-current  curves 
versus  flow-field  pressure  drop  of  a  PEM  fuel  cell  by  experimen¬ 
tal  study.  Grujicic  et  al.  [11]  combined  a  steady-state,  single-phase 
electrochemical  model  with  a  non-linear  constrained  optimization 
procedure  to  maximize  the  performance  of  the  cathode  in  a  PEMFC. 
Chen  et  al.  [12]  presented  a  study  of  the  internal  and  external  struc¬ 
ture  effects  on  the  fuel  cell  steady  and  transient  operation.  Cheng 
et  al.  [13]  integrated  a  simplified  conjugate-gradient  method  with 
CFD  solver  to  determine  the  optimum  set  of  geometric  parameters. 
Mo  et  al.  [14]  reported  an  optimization  study  of  PEMFCs  using  a 
hybrid  genetic  algorithm. 


The  major  challenge  to  geometry  optimization  in  computational 
fuel  cell  dynamics-based  design  lies  on  generating  a  high  qual¬ 
ity  computational  mesh  and  subsequently  running  the  simulation 
to  obtain  the  objective  function  through  automated  evaluation 
with  commercial  software.  The  COMSOL™  is  well  integrated  with 
MATLAB™  and  can  be  implemented  in  MATLAB™  as  a  toolbox, 
which  offers  a  useful  option  for  the  self-automated  geometry  opti¬ 
mization  process  in  PEMFCs.  This  study  focuses  on  the  geometry 
optimization,  in  which  the  effects  of  the  anode  channel  ratio  and  the 
cathode  open  ratio  on  the  performance  of  the  air-breathing  PEMFC 
are  investigated  in  order  to  improve  the  cell  performance.  The 
sequential  quadratic  programming  (SQP)  method  [15-17],  being  a 
high  efficiency  method,  is  selected  to  deal  with  the  constrained 
design  problem,  while  the  objective  function  is  calculated  by  run¬ 
ning  the  3D  simulation  script  of  COMSOL™  under  the  MATLAB™ 
environment. 


2.  Modeling  with  COMSOL™ 

Modeling  of  an  air-breathing  PEMFC  is  a  challenge  since  the  pro¬ 
cess  is  a  multi-physics  phenomenon  that  involves  electrochemical 
reactions,  flow  dynamics  in  the  gas  channels,  diffusion  and  convec¬ 
tion  in  the  GDL  layer,  and  heat  and  mass  transportation  in  porous 
media.  In  this  study,  a  steady-state,  single-phase  three-dimensional 
model  for  an  air-breathing  PEMFC  is  developed  using  COMSOL™ 
to  simulate  the  performance  of  the  cell  which  is  considered  as  the 
objective  function  value  for  the  optimization  process. 

Fig.  1(a)  schematically  shows  an  air-breathing  PEMFC  and  its 
various  components  including  the  membrane,  the  flow  channels, 
GDLs  and  catalyst  layers  on  both  the  anode  and  the  cathode  sides. 
The  computational  geometry  is  shown  in  Fig.  1(b).  The  catalyst 
layers  at  both  the  anode  and  cathode  are  simplified  as  bound¬ 
aries  between  the  GDLs  and  the  membrane,  with  the  thicknesses 
of  the  catalyst  layers  accounted  in  the  membrane’s  dimension. 
The  numerical  domain  computed  is  based  on  a  repetitive  unit  of 
the  PEMFC  between  the  channel  and  the  rib,  and  thus  the  exter¬ 
nal  boundaries  of  the  domain  can  be  treated  as  symmetric.  Air  is 
induced  into  the  cathode  channels  by  virtue  of  natural  convection, 
whereas  hydrogen  is  delivered  into  the  anode  channels  through 
forced  convection.  The  pressures  of  both  the  air  and  hydrogen  are 
set  at  one  standard  atmosphere. 

Here,  the  applicability  of  the  continuum  theory  needs  to  be 
checked  since  the  characteristic  size  of  the  computational  domain  is 
very  small.  It  is  well  known  that  the  Knudsen  number  (Kn),  defined 
as  Kn  =  A/L,  is  used  to  measure  the  ratio  of  the  molecular  mean 
free  path  length  (A)  to  a  representative  physical  length  scale  (L). 
The  condition  Kn  >  1  is  called  the  free  molecule  regime,  and  for 
Kn  <$c  1  continuum  theory  applies.  For  the  particular  simulation  in 
this  study,  the  continuum  theory  is  valid  since  the  Knudsen  number 
( Kn )  is  at  the  order  of  10-5.  Assumptions  used  in  the  model  include: 


(1 )  the  Water  exists  in  the  fuel  cell  as  vapour 

(2)  the  electrodes  and  membrane  are  made  of  homogeneous  mate¬ 
rials  with  uniform  morphological  properties 

(3)  the  cell  is  considered  to  be  isothermal,  i.e.,  the  temperature 
distribution  across  the  cell  is  uniform 

(4)  humid  air  is  considered  on  the  cathode  side  while  only  hydro¬ 
gen  and  water  vapour  are  considered  on  the  anode  side. 


Based  on  the  model  assumptions  and  the  scale  analyses,  the 
flow-field  of  reactant  gases  in  the  gas  channel  is  governed  by 
the  continuity  equation  and  ensures  the  mass  conservation  and 
the  steady  state  incompressible  Navier-Stokes  equation  which 
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Fig.  1.  Scheme  of  (a)  air-breathing  PEMFC  fuel  cell  and  (b)  computational  domain. 


describes  the  momentum  conservation: 


V  •  u  =  0  (1) 

pS=-vp+F  (2) 

where  u  is  the  velocity  vector,  p  stands  for  pressure,  F  is  the  body 
force,  and  p  is  the  mixture  density. 

The  multi-component  diffusion  and  convection  phenomena  are 
described  by  the  Maxwell-Stafen  equation,  which  has  the  general 
form: 
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where  Dy  is  the  diffusion  coefficient,  DT  is  the  multi-component 
thermal  diffusion  coefficient,  Q  is  the  source  term,  x  is  the  species 


mole  fraction,  go  is  the  species  mass  fraction,  T  is  the  tempera¬ 
ture,  and  M  is  the  molecular  mass.  Temperature-driven  diffusion  is 
insignificant  due  to  an  assumed  uniform  temperature  and  therefore 
the  source  term  Q*  is  set  to  zero. 

On  the  cathode  side,  two  transport  equations  are  solved  to  attain 
the  mass  fractions  of  oxygen  and  water.  The  mass  fraction  of  nitro¬ 
gen  can  be  obtained  from  the  mass  balance  equation: 

&>n2  =  1  o)o2  -  coH2o  (4) 

Similarly,  on  the  anode  side,  the  transport  equation  of  hydrogen 
can  be  solved  and  the  mass  fraction  of  water  obtained  by 

&>h2o  =  1  -  &>h2  (5) 

Since  the  GDLs  are  porous  media,  the  velocity  distribution  is 
therefore  formulated  by  Darcy’s  law: 

u  —  -^Vp  (6) 

where  /c  is  the  permeability  and  p  is  the  dynamic  viscosity. 

At  the  electrode-membrane  boundary,  the  mass  fluxes  of  hydro¬ 
gen  in  the  anode  (Nh2  ),  and  of  oxygen  (Nq2  )  and  water  (Nh2o  )  in  the 


Fig.  2.  Computational  domain  (a)  top  view  and  (b)  3D  view. 


Table  1 

Geometry  and  functional  parameters. 
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Description 

Value 

Anode  GDL  thickness 

0.35  mm 

Cathode  GDL  thickness 

0.35  mm 

Thickness  of  the  membrane  and  catalyst  layers 

0.035  mm 

Height  of  the  anode  channel 

0.8  mm 

Height  of  the  cathode  channel 

1.5  mm 

Anode  exchange  current  density 

1000  Am-2 

Cathode  exchange  current  density 

0.1  Am-2 

Faraday’s  constant 

96.487C  mol-1 

Membrane  conductivity 

0.9  S  cm-1 

Solid-phase  conductivity 

1250  8  cm-1 

Porosity  of  GDL 

0.35 

Porosity  of  catalyst  layer 

0.1 

cathode,  are  determined  by  the  following  reaction  rates,  respec¬ 
tively: 


Nh2  =  -^MH2,  No2=-Mmo2,  Nh2o^m„20  (7) 

where  ja  and  jc  are  the  transfer  current  densities  corresponding 
to  the  electrochemical  reaction  at  the  anode  and  cathode  catalyst 
layers. 

The  agglomerate  model  describes  the  current  density  in  an 
active  layer  consisting  of  agglomerates  of  ionic  conductor  material 
and  electrically  conducting  particles  covered  partially  with  cata¬ 
lyst.  Governing  equations  for  the  current  density  in  the  anode  and 
cathode  are 
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where  Dagg  is  the  agglomerate  gas  diffusivity,  1?  is  the  gas  con¬ 
stant,  ftagg  is  the  agglomerate  radius,  s  is  the  specific  area  of  the 
catalyst  inside  the  agglomerate,  F  is  the  Faraday’s  constant,  cref 
and  cagg  are  the  reference  concentrations  of  the  species  and  that 
in  the  agglomerate  surface,  while  J0  is  the  exchange  current  den- 
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Table  2 

Geometry  and  functional  parameters. 


Description 

Value 

Anode  GDL  thickness 

0.25  mm 

Cathode  GDL  thickness 

0.25  mm 

Thickness  of  the  membrane  and  catalyst  layers 

0.25  mm 

Height  of  the  anode  channel 

0.5  mm 

Height  of  the  cathode  channel 

1.2  mm 

Anode  exchange  current  density 

1000  Am-2 

Cathode  exchange  current  density 

0.1  Am-2 

Faraday’s  constant 

96,487  C  mol-1 

Membrane  conductivity 

9  8  cm-1 

Solid-phase  conductivity 

1000  S  cm-1 

Permeability 

10“13  m2 

Anode  electrode  potential 

ov 

sity,  and  T  is  the  temperature.  Finally,  the  overvoltage  rjVt  is  defined 
as  rj\j  =  0s  -0m  -£eq,  where  0S  and  0m  are  the  potentials  in  the 
electrode  and  in  the  membrane,  respectively.  Eeq  is  the  equilibrium 
voltage. 

The  conductive  media  direct  current  application  mode  in 
COMSOL™  describes  the  potential  distributions  in  the  solid-phase 
GDL  and  the  membrane  using  the  following  equations: 

v  •  (-tfs.eff^s)  =  Ss,  V  •  (-(Xm,effV0m)  =  Sm  (9) 

where  crs  eff  is  the  solid-phase  effective  electronic  conductivity,  and 
crm,eff  is  the  membrane  ionic  conductivity.  5S  and  Sm  are  the  current 
source  terms  in  the  solid  phase  and  the  electrolyte  phase,  respec¬ 
tively. 


Fig.  3.  Polarization  curves. 


Fig.  5.  Optimization  convergence  histories. 
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3.  Model  validation 

Due  to  the  limitation  on  the  extensive  memory  requirement  of 
COMSOL™,  only  a  small  section  shown  in  Fig.  2(a)  is  selected  to 
simulate  the  behaviour  of  the  PEMFC.  Thus,  the  computation  is  lim¬ 
ited  to  a  repetitive  model  between  the  channel  and  the  rib.  The 
external  boundaries  can  be  treated  as  symmetric.  The  3D  simula¬ 
tion  domain  is  schematically  outlined  in  Fig.  2(b),  which  consists 
of  the  cathode  air-breathing  slot,  two  GDLs,  the  membrane  and  the 
catalyst  layers,  and  the  anode  channel  for  the  hydrogen  supply.  The 
cathode  open  ratio  is  selected  to  be  80%.  In  order  to  make  sure 
that  the  boundaries  can  be  treated  as  symmetric,  the  widths  of  the 
channel  and  the  rib  are  selected  as  1.6  and  0.4  mm  (4:1),  respec¬ 
tively,  at  the  cathode  side,  and  1  and  1  mm  (1 :1 )  for  the  anode  side, 
instead  of  4  and  1  mm  for  the  cathode  as  designed  by  Schmitz  et  al. 
[18],  whose  experimental  results  are  selected  to  validate  the  cur¬ 


rent  model.  All  other  geometric  parameters  such  as  the  thicknesses 
of  the  GDLs  and  membrane,  the  depths  of  the  air-breathing  chan¬ 
nel  and  the  hydrogen  anode  channel,  and  so  on,  are  represented 
by  those  described  by  Schmitz  et  al.  [18]  and  are  listed  in  Table  1. 
Fuel  cell  temperatures  vary  under  different  current  conditions  and 
Schmitz  et  al.  [18]  only  reported  those  corresponding  to  two  con¬ 
ditions,  namely,  63  °C  for  the  cell  under  short-circuit  conditions 
and  43  °C  at  the  maximum  power  point  with  the  opening  ratio  of 
80%.  The  temperature  in  this  study  is  selected  as  43  °C  according 
to  the  maximum  power  point  condition  in  Schmitz  et  al.  [18].  An 
unstructured  mesh  with  3160  elements  is  used  for  this  simulation. 
Fig.  3  shows  the  polarization  curves  obtained  from  the  simulation 
and  that  from  the  experimental  study  by  Schmitz  et  al.  [18]  at  a 
cathode  open  ratio  of  80%.  It  is  found  that  the  simulated  cell  perfor¬ 
mance  at  the  maximum  power  point  (43  °C,  0.4  V)  is  in  reasonable 
agreement  with  the  experimental  results  of  Schmitz  et  al.  [18], 


Min:  0.177 


(b)  Mass  fraction  w02  Max:  0.230 


Min:  0.177 


Min:  0.177 


Fig.  6.  Mass  fraction  distributions  of  oxygen  with  open  ratios  of  (a)  15%,  (b)  49.8%  and  (c)  80%. 
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indicating  that  the  model  proposed  here  can  generate  relatively 
accurate  data. 

4.  Grid  size  independency  study 

The  sensitivity  of  the  converged  fuel  cell  simulation  to  the  size 
of  the  computational  mesh  is  demonstrated  by  carrying  out  the 
simulation  under  the  same  operating  conditions  and  consistent 
functional  parameters  for  an  air-breathing  PEMFC,  but  with  dif¬ 
ferent  mesh  sizes.  The  grid  size  independence  studies  have  been 
carried  out  with  three  different  mesh  sizes  with  3085,  6811  and 
11,764  elements,  respectively.  Fig.  4  shows  the  current  densities 
obtained  for  the  three  different  mesh  sizes  based  on  a  cell  operat¬ 
ing  at  a  316  K,  a  cell  voltage  of  0.7  V,  and  with  the  cathode  open  ratio 
set  at  80%  (33.2136,  33.2204,  and  33.2244  mA  cm-2  corresponding 
to  a  mesh  size  of  3085,  6811  and  11,764,  respectively).  The  current 


densities  obtained  with  different  element  numbers  exhibit  little 
difference  and  converge  to  a  limiting  value  as  the  number  of  mesh 
element  increases.  This  suggests  that  using  the  coarse  mesh  with 
3085  elements  not  only  provides  reasonably  accurate  results,  but 
also  saves  computational  costs  significantly,  especially  when  a  con¬ 
siderable  number  of  simulation  iterations  are  required  to  evaluate 
the  objective  function  values  needed  for  the  optimization  pro¬ 
cesses.  Therefore,  the  same  coarse  size  mesh  is  used  to  evaluate  the 
objective  function  values  for  the  geometry  optimization  studies. 

5.  Optimization  method 

Sequential  quadratic  programming  methods  (SQP)  have  been 
developed  to  solve  reliably  constrained  optimization  problems 
with  multiple  variables  and  constraints.  These  methods  require 
only  a  few  evaluations  of  the  objective  functions  and  have  shown 
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Fig.  7.  Mass  fraction  distributions  of  hydrogen  with  open  ratios  of  (a)  15%,  (b)  49.8%  and  (c)  80%. 
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their  ability  in  terms  of  accuracy  and  percentage  of  successful  solu¬ 
tions. 

The  general  method  is  stated  here  briefly.  Given  a  problem 
described  as 


minmiz  e/(x) 
subject  to 

&z(x)  =  0  i  =  1 , . . . ,  me 
gj(x)  <0  i  =  me  +  1 , . . . ,  m 


(10) 


where,  x  is  the  p-dimensional  vector  of  the  design  variables, /(x)  is 
the  objective  function. 

The  basic  idea  of  the  SQP  method  is  to  formulate  a  quadratic  pro¬ 
gramming  (QP)  sub-problem  based  on  a  quadratic  approximation 
of  the  Lagrangian  function: 

m 

L(x,X)=f(x)  +  yfxigi(x)  (11) 

1=1 

where  Lagrange  multipliers  Xif  i=  1,  m  are  used  to  balance  the 
deviations  in  magnitude  of  the  objective  function  and  constraint 
gradients.  Only  active  constraints  are  included  in  this  operation, 
constraints  that  are  not  active  excluded.  Hence,  Eq.  (10)  can  be  sim¬ 
plified  by  assuming  that  bound  constraints  have  been  expressed  as 


inequality  constraints,  and  the  QP  sub-problem  is  formed  as 

minmize^drHfcd  +  V/(xfc)Td 

d  e  mn  2 

subject  to  (i2) 

Vg,(xk  fd  +  gj(xk)  =  0  i  =  1, ...  ,me 
Vgj(xk)Td  +  gj(xk)  <0  i  =  me  +  1, ...  ,m 

where  H  is  a  positive  definite  approximation  of  the  Hessian  matrix 
of  the  Lagrangian  function  (Eq.  (11 )),  and  can  be  updated  by  any  of 
the  quasi-Newton  methods,  d  is  the  vector  of  the  design  variables 
of  the  sub-problem. 

The  solution  of  the  sub-problem  is  used  to  form  a  new  iteration: 

Xfc+1  =xk+akdk  (13) 

The  step  size  parameter  ak  is  determined  by  an  appropriate  line 
search  procedure  so  that  a  sufficient  decrease  in  a  merit  function  is 
obtained.  dk  stands  for  the  search  direction. 

The  SQP  implementation  consists  of  three  main  steps  namely: 
updating  the  Hessian  matrix,  solving  the  quadratic  programming 
sub-problem,  and  forming  the  new  iteration.  The  implementation 
details  will  not  be  discussed  here  but  can  be  found  in  [15-17]. 

In  this  paper,  the  optimization  problem  is  solved  using  the  SQP 
algorithm  in  the  MATLAB™  optimization  toolbox  to  minimize  the 


Min:  0.0306  Min:  0.0306 


Max:  0.230 


Min:  0.0306 


Fig.  8.  Mass  fraction  distributions  of  oxygen  with  open  ratios  of  (a)  15%,  (b)  49.5%  and  (c)  80%. 
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(a)  Slice:  Mass  fraction,  wH2 
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Fig.  9.  Mass  fraction  distributions  of  hydrogen  with  open  ratios  of  (a)  15%,  (b)  49.5%  and  (c)  80%. 


Fig.  10.  Polarization  curves  with  three  different  open  ratios. 


Fig.  11.  Optimization  convergence  history. 
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Fig.  12.  Mass  fraction  distributions  of  oxygen  with  channel  ratios  of  (a)  20%,  (b)  34.7%  and  (c)  80%. 


negative  of  the  average  current  density  at  the  cathode  side.  The 
objective  function  is  evaluated  by  running  the  3D  COMSOL™  sim¬ 
ulation  script  under  the  MATLAB™  environment. 

6.  Results  and  discussion 

Two  cases  are  studied  to  find  out  the  optimum  geometric  design 
of  the  air-breathing  fuel  cell: 

Case  1.  Keeping  the  anode  channel  size  constant  while  optimizing 
the  cathode  open  ratio. 

Case  2.  Fixing  the  cathode  open  ratio  while  locating  the  optimum 
anode  channel  ratio. 

For  the  above  two  cases,  the  fuel  cell  operating  temperature 
is  maintained  at  353  K  and  the  pressure  is  kept  at  one  standard 
atmosphere.  Cell  design  parameters,  material  properties  and  the 
functional  parameters  used  in  the  simulations  are  listed  in  Table  2. 
The  objective  function  is  set  as  the  current  density,  while  the  design 
variables  are  selected  as  the  cathode  open  ratio  for  Case  1  and  the 
anode  channel  ratio  for  Case  2,  respectively.  Pure  hydrogen  is  sup¬ 
plied  at  the  anode  side  and  humidified  air  is  supplied  to  the  cathode 
side. 


6.1.  Case  1:  optimization  of  cathode  open  ratio 

Primary  studies  indicate  that  a  small  anode  channel  ratio  can 
supply  sufficient  hydrogen  for  a  cell  with  a  relatively  large  cathode 
open  ratio  when  the  cell  operates  at  the  temperature,  pressure,  and 
the  functional  parameters  listed  in  Table  2.  Therefore,  in  Case  1,  the 
anode  channel  ratio  is  fixed  at  10%,  and  the  cathode  open  ratio  is  set 
as  the  independent  design  variable.  Lower  and  upper  boundaries  of 
the  design  variable  are  set  at  10%  and  90%,  respectively.  The  current 
density  is  considered  as  the  objective  function.  The  values  of  the 
operating  cell  voltages  of  the  fuel  cell  are  selected  as  0.7  and  0.4  V  in 
the  optimization  process  to  determine  the  corresponding  optimum 
cathode  open  ratios,  respectively. 

Since  the  SQP  approach  is  a  local  optimization  method,  there 
are  possibilities  that  the  optimization  process  converges  to  a  local 
optimum  instead  of  a  global  result.  In  order  to  make  sure  that  the 
optimization  converges  to  a  global  optimal  solution,  the  searching 
space  is  divided  into  a  few  sub-domains  to  guarantee  that  the  whole 
spectrum  of  design  variables  (10-90%)  has  been  searched. 

Convergence  histories  with  cell  voltages  at  0.7  and  0.4  V  are 
presented  in  Fig.  5  with  square  symbols  and  filled  delta  symbols, 
respectively.  The  current  density  increases  with  enlarging  the  open 
ratio  from  10%  to  49.8%,  and  beyond  49.8%  there  is  no  significant 


X.Q..  Xing  et  al.  /  Journal  of  Power  Sources  186  (2009)  10-21 


19 


increment  in  the  current  density.  Therefore,  the  optimum  current 
density  occurs  at  the  open  ratio  of  49.8%  at  a  cell  voltage  of  0.7  V. 
The  optimum  open  ratio  is  found  to  be  49.5%  when  the  cell  voltage 
is  fixed  at  0.4  V.  The  convergence  histories  show  the  same  trends 
for  the  two  different  cell  voltages  considered,  i.e.,  the  current  den¬ 
sity  rises  with  the  open  ratio  increasing  from  10%  to  the  optimum 
open  ratio  value.  No  obvious  increment  in  the  current  density  is 
observed  with  further  increases  in  the  open  ratio  beyond  the  opti¬ 
mum. 

For  a  cell  voltage  of  0.7  V,  mass  fraction  distributions  of  the  oxy¬ 
gen  at  the  cathode  side  and  the  hydrogen  at  the  anode  side  are 
presented  in  Figs.  6  and  7  with  open  ratios  of  15%,  49.8%,  and  80%, 
respectively.  From  Fig.  6,  it  can  be  clearly  seen  that  oxygen  increases 
when  the  cathode  open  ratio  rises  from  15%  to  49.8%,  and  further  to 
80%  since  more  air  is  supplied  to  the  cell.  The  data  in  Fig.  7  indicates 
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that  when  the  anode  channel  ratio  is  fixed  at  10%,  the  hydrogen 
concentration  is  slightly  affected  by  the  cathode  open  ratio  vary¬ 
ing  from  15%  to  80%.  This  is  due  to  the  fact  that  as  the  cathode 
open  ratio  is  at  15%,  the  oxygen  supply  would  become  the  limit¬ 
ing  factor  to  support  the  chemical  reaction,  leading  to  over-supply 
of  hydrogen.  Flence,  the  mass  fraction  of  hydrogen  is  seen  to  be 
slightly  higher  than  that  of  the  other  two  open  ratio  cases  while  the 
mass  fraction  of  the  oxygen  at  the  cathode  is  seen  to  be  the  low¬ 
est.  As  the  cathode  open  ratio  is  increased,  oxygen  supply  becomes 
more  abundant,  resulting  in  a  higher  consumption  of  hydrogen, 
and  leading  to  a  higher  current  density.  Further  increase  in  the 
open  ratio  beyond  49.8%  results  in  no  obvious  improvement  in  cell 
performance,  which  suggests  that  when  the  open  ratio  is  beyond 
49.8%,  hydrogen  supply  and  oxygen  supply  are  both  sufficient,  and 
the  chemical  reaction  reaches  a  balance  under  the  given  cell  volt- 
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Fig.  13.  Mass  fraction  distributions  of  hydrogen  with  channel  ratios  of  (a)  20%,  (b)  34.7%  and  (c)  80%. 
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age.  Consequently,  no  further  increment  in  the  current  density  is 
obtained. 

For  a  cell  voltage  of  0.4  V,  mass  fraction  distributions  of  the 
oxygen  at  the  cathode  side  and  hydrogen  at  the  anode  side  are  pre¬ 
sented  in  Figs.  8  and  9  with  open  ratios  of  15%,  49.5%,  and  80%, 
respectively.  Qualitatively,  a  similar  trend  to  that  observed  at  0.7  V 
is  found  for  both  oxygen  and  hydrogen  distributions  in  consider¬ 
ation  at  the  various  cathode  open  ratios.  Flowever,  quantitatively, 
the  0.4  and  0.7  V  cases  exhibit  very  different  phenomena,  as  seen 
by  comparing  the  mass  fraction  values  of  the  hydrogen  and  the 
oxygen  left  in  the  GDLs  and  the  channels  or  the  air-breathing  slots. 
For  a  cell  voltage  of  0.7  V,  the  mass  fractions  of  the  oxygen  and  the 
hydrogen  left  in  the  GDLs  and  the  air-breathing  slots  or  the  chan¬ 
nels  are  quite  high  when  the  cathode  open  ratio  is  larger  or  equal 
to  49.8%,  as  shown  in  Figs.  6  and  7.  Flence  for  this  case,  cell  per¬ 
formance  is  limited  by  the  rate  of  the  chemical  reaction.  For  a  cell 
voltage  of  0.4  V,  the  data  in  Figs.  8  and  9  show  that  the  mass  frac¬ 
tions  of  hydrogen  left  in  the  GDL  and  the  channel  at  anode  are  very 
low  when  the  open  ratio  is  larger  or  equal  to  49.5%.  This  suggests 
that  no  further  improvement  in  cell  performance  is  obtained  when 
the  open  ratio  is  beyond  49.5%  since  the  hydrogen  is  depleted  by 
the  chemical  reaction.  Hence,  for  this  case,  the  cell  performance  is 
limited  by  the  hydrogen  supply. 

The  polarization  curves  with  different  open  ratios  of  15%,  49.8% 
and  80%  are  shown  in  Fig.  10.  It  is  clearly  noticeable  that  the  cell 
performance  is  very  low  when  the  open  ratio  is  as  small  as  15%  due 
to  insufficient  oxygen  supply.  The  difference  in  cell  performance 
between  49.8%  and  80%  is  not  significant  at  high  operating  voltages 
(i.e.,  larger  or  equal  to  0.4  V).  When  the  operating  voltages  are  lower 
than  0.4  V,  with  an  80%  cathode  open  ratio,  the  hydrogen  supply  is 
far  from  sufficient  with  an  anode  channel  ratio  of  10%,  and  conse¬ 
quently,  maximum  current  occurs  at  370.4  mA  cm-2,  which  is  even 
lower  than  the  maximum  when  the  cathode  open  ratio  is  49.8%. 

6.2.  Case  2:  optimization  of  channel  ratio  at  anode 

In  this  case,  the  anode  channel  ratio  is  designed  as  the  inde¬ 
pendent  variable  and  the  current  density  is  still  considered  as  the 
objective  function.  The  operating  cell  voltage  of  the  fuel  cell  is 
selected  as  0.7  V  in  the  optimization  process  to  determine  the  cor¬ 
responding  optimum  anode  channel  ratio. 

For  the  same  reason  previously  mentioned  in  Case  1,  a  relatively 
small  anode  channel  ratio  can  supply  sufficient  hydrogen  for  a  cell 
with  a  large  cathode  open  ratio,  especially  when  the  cell  voltage  is 
as  high  as  0.7  V.  Hence,  in  Case  2,  the  cathode  open  ratio  is  fixed 


Fig.  14.  Polarization  curves  with  three  different  channel  ratios. 


at  80%.  Lower  and  upper  boundaries  of  the  design  variable  (anode 
channel  ratio)  are  set  at  20%  and  80%,  respectively.  The  convergence 
history  of  the  optimization  process  is  presented  in  Fig.  11.  The  opti¬ 
mization  history  indicates  that  the  optimum  anode  channel  ratio 
is  located  at  34.7%.  The  variations  of  the  current  density  obtained 
with  different  anode  channel  ratios  are  very  small.  Mass  fraction 
distributions  of  the  oxygen  in  the  GDLs  and  the  air-breathing  slots 
are  presented  in  Fig.  12.  The  mass  fraction  of  the  oxygen  is  seen  to 
differ  insignificantly  at  the  different  anode  channel  ratios  of  20%, 
34.7%  and  80%,  and  the  minimum  mass  fraction  of  the  oxygen  left  in 
the  GDLs  and  the  slots  is  0.218,  which  is  94.8%  of  the  mass  fraction 
of  the  oxygen  in  the  ambient.  Fig.  13  shows  the  mass  fraction  distri¬ 
butions  of  hydrogen  at  the  anode  side.  It  can  be  seen  that  hydrogen 
is  present  in  a  high  amount  even  when  the  channel  ratio  is  as  low 
as  20%.  The  optimization  history,  together  with  the  mass  fraction 
distributions  of  the  oxygen  at  the  cathode  side  and  the  hydrogen 
at  the  anode  side,  indicate  that  the  fuel  supply  is  in  abundance. 
The  chemical  reaction  rate  changes  little  with  the  channel  ratios, 
and  therefore  results  in  a  flat  curve  presented  in  the  optimization 
history. 

Fig.  14  presents  the  three  polarization  curves,  corresponding  to 
the  anode  channel  ratios  of  20%,  34.7%,  and  80%,  respectively.  There 
is  insignificant  difference  among  the  three  polarization  curves 
when  the  cell  voltage  is  larger  than  0.2  V.  Below  0.2  V,  the  hydro¬ 
gen  supply  is  insufficient  with  an  anode  channel  ratio  of  20%.  The 
higher  the  anode  channel  ratio,  the  higher  the  current  density. 

Studies  of  the  mass  fraction  distributions  indicate  that  the 
supplies  of  hydrogen  and  oxygen  are  sufficient  for  the  chemical 
reaction  when  the  cathode  open  ratio  is  fixed  at  80%,  the  anode 
channel  ratio  is  greater  than  20%,  and  the  cell  voltage  is  larger  than 
0.2  V.  The  effect  of  the  anode  channel  ratios  on  the  chemical  reac¬ 
tion  rate  is  negligible,  leading  to  no  significant  difference  among 
the  three  polarization  curves. 

7.  Conclusion 

Using  the  deep  integration  between  COMSOL™  and  MATLAB™, 
an  automated  geometry  optimization  of  an  air-breathing  PEMFC, 
searching  for  the  optimum  anode  channel  ratio  and  the  cathode 
open  ratio  for  the  enhanced  fuel  cell  performance,  is  investigated. 
Simulation  results  for  the  air-breathing  PEMFC  operating  at  353  K 
and  one  standard  atmosphere  pressure  show  the  following: 

(1)  When  the  anode  open  ratio  is  fixed  at  10%,  the  optimum  cath¬ 
ode  open  ratios  at  two  different  operation  voltages  of  0.7  V  and 
0.4  V  are  rather  similar,  namely,  49.8%  with  a  cell  voltage  of 
0.7  V  and  49.5%  with  a  cell  voltage  of  0.4  V.  Polarization  curves 
indicate  that  when  the  cathode  open  ratio  is  less  than  15%,  the 
oxygen  supply  is  insufficient,  and  results  in  a  low  cell  perfor¬ 
mance.  When  the  operating  voltage  is  larger  than  or  equal  to 
0.4  V,  there  is  insignificant  difference  in  the  cell  performance 
between  the  cathode  open  ratios  of  49.8%  and  80%.  When  oper¬ 
ation  voltages  are  lower  than  0.4  V  (near  to  the  limiting  current 
density),  however,  with  an  80%  cathode  open  ratio,  the  hydro¬ 
gen  supply  will  be  depleted  with  a  small  anode  channel  ratio  of 
10%,  and  consequently,  the  cell  performs  worse  than  that  with 
a  cathode  open  ratio  of  49.8%. 

(2)  When  the  cathode  open  ratio  is  set  at  80%,  the  optimum  anode 
channel  ratio  is  found  to  be  34.7%  at  a  cell  voltage  of  0.7  V.  The 
optimization  history,  together  with  the  mass  fraction  distribu¬ 
tions  of  the  oxygen  at  the  cathode  side  and  the  mass  fraction 
distributions  of  the  hydrogen  at  the  anode  side,  suggest  that 
the  fuel  supply  is  in  abundance.  The  rate  of  chemical  reaction 
varies  insignificantly  with  the  anode  channel  ratio,  resulting  in 
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an  asymptotic  curve  presented  in  the  optimization  history  and 
little  difference  among  the  three  polarization  curves. 
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